!>\brief
!! In this module the local matrix and right hand side for the CG
!! method are computed.
!!
!! \n
!!
!! The local matrix does not include any static condensation; it
!! includes the Neumann type boundary conditions but not the Dirichlet
!! type ones.
!!
!! \note There are no consistency checks on the coefficients returned
!! by the \c coeff* functions of \c mod_testcases, which thus have to
!! be correct.
!<----------------------------------------------------------------------
module mod_cg_locmat

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_master_el, only: &
   mod_master_el_initialized, &
   t_me

 use mod_base, only: &
   mod_base_initialized, &
   t_base

 use mod_grid, only: &
   mod_grid_initialized, &
   t_e, affmap

 use mod_bcs, only: &
   mod_bcs_initialized, &
   p_t_b_v, p_t_b_s, p_t_b_e, &
   b_dir,   b_neu,   b_ddc
 
 use mod_testcases, only: &
   mod_testcases_initialized, &
   coeff_diff,  &
   coeff_adv,   &
   coeff_re,    &
   coeff_f,     &
   coeff_neu,   &
   coeff_rob

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_cg_locmat_constructor, &
   mod_cg_locmat_destructor,  &
   mod_cg_locmat_initialized, &
   cg_locmat, &
   t_bcs_error

 private

!-----------------------------------------------------------------------

! Module types and parameters

 ! This type is used to indicate that an error has occurred in the
 ! computation of the boundary condition
 type t_bcs_error
   logical :: lerr
   character(len=1000) :: message
 end type t_bcs_error

! Module variables

 ! public members
 logical, protected ::               &
   mod_cg_locmat_initialized = .false.
 ! private members
 character(len=*), parameter :: &
   this_mod_name = 'mod_cg_locmat'

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_cg_locmat_constructor()

  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized.eqv..false.) .or. &
       (mod_kinds_initialized.eqv..false.)    .or. &
       (mod_master_el_initialized.eqv..false.).or. &
       (mod_base_initialized.eqv..false.)     .or. &
       (mod_grid_initialized.eqv..false.)     .or. &
       (mod_bcs_initialized.eqv..false.)      .or. &
       (mod_testcases_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_cg_locmat_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   mod_cg_locmat_initialized = .true.
 end subroutine mod_cg_locmat_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_cg_locmat_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_cg_locmat_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_cg_locmat_initialized = .false.
 end subroutine mod_cg_locmat_destructor

!-----------------------------------------------------------------------
 
 !> Computation of the local matrix and right hand side.
 !!
 !! The local matrix is
 !! \f{displaymath}{
 !!  A^K_{ij} = \int_K \left[(\mu\nabla\varphi_j - a\varphi_j) \cdot
 !!  \nabla\varphi_i + \sigma\varphi_j\varphi_j\right]\,dx
 !! \f}
 !! and the right hand side is
 !! \f{displaymath}{
 !!  f^K_{i} = \int_K f\varphi_i\,dx + bcs.
 !! \f}
 !<
 pure subroutine cg_locmat(aak,fk,base,e,be,bcs_err)
  real(wp), intent(out) :: aak(:,:), fk(:)
  type(t_base), intent(in) :: base
  type(t_e), intent(in) :: e
  type(p_t_b_e), intent(in) :: be
  type(t_bcs_error), intent(out) :: bcs_err
 
  integer :: l, i, j
  real(wp) :: &
    mu(e%m,e%m,base%m), a(e%m,base%m),        &
    sigma(base%m), f(base%m), bi_mu_bit(e%d,e%d), &
    bi_a(e%d), xg(e%m,base%m), mu_grad(e%m), &
    a_grad_sigma_p
  character(len=*), parameter :: &
    this_sub_name = 'cg_locmat'

   ! Gaussian nodes in physical space
   xg = affmap(e,base%xig)

   ! problem coefficients
   mu    = coeff_diff(xg)
   a     = coeff_adv(xg)
   sigma = coeff_re(xg)
   f     = coeff_f(xg)

   ! local matrix including diffusion, advection and reaction
   aak = 0.0_wp
   fk  = 0.0_wp
   do l=1,base%m
     ! bi_mu_bit = b^{-1} * mu * b^{-T}
     bi_mu_bit = matmul( e%bi , matmul(mu(:,:,l),transpose(e%bi)) )
     ! bi_a = b^{-1} * a
     bi_a = matmul( e%bi , a(:,l) )
     do i=1,base%pk
       mu_grad = matmul( bi_mu_bit , base%gradp(:,i,l) )
       a_grad_sigma_p = - dot_product( bi_a , base%gradp(:,i,l) ) &
                                  + sigma(l) * base%p(i,l)
       do j=1,base%pk
         aak(i,j) = aak(i,j) + base%wg(l) * (                        &
                   dot_product( base%gradp(:,j,l) , mu_grad )        &
                                    + base%p(j,l) * a_grad_sigma_p )
       enddo
       fk(i) = fk(i) + base%wg(l) * f(l) * base%p(i,l)
     enddo
   enddo
   aak = e%det_b*aak
   fk  = e%det_b*fk

   ! add boundary fluxes if present
   bcs_err%lerr = .false.
   if(associated(be%p)) then
     call cg_boundary_flux(aak,fk,base,be,bcs_err)
   endif
   !--------------------------------------------------------------------
 
 end subroutine cg_locmat
 
!-----------------------------------------------------------------------
 
 pure subroutine cg_boundary_flux(aak,fk,base,be,err)
  real(wp), intent(inout) :: aak(:,:), fk(:)
  type(t_base), intent(in) :: base
  type(p_t_b_e), intent(in) :: be
  type(t_bcs_error), intent(out) :: err

  integer :: isb, isl, breg, l, i, j
                              !integer l_e(base%me)
  real(wp) :: xgs(be%p%e%m,base%ms), a(be%p%e%m,base%ms), &
    gamma(base%ms), hb(base%ms), an(base%ms), ap(base%ms),&
    wgsa(base%ms)
 
   err%lerr = .false.
   do isb=1,be%p%nbs ! loop on the boundary sides
     btype: select case(be%p%bs(isb)%p%bc)
      case(b_dir,b_ddc)
       ! nothing to do in this subroutine
      case(b_neu)
 
       isl = be%p%bs(isb)%p%s%isl(1) ! local index of the side
       ! Gaussian nodes in physical space
       xgs = affmap(be%p%e,base%xigb(:,:,isl))
       ! side weigths reordered and scaled by the side area
       wgsa = (be%p%e%s(isl)%p%a/base%me%voldm1) * &
                   base%wgs(base%stab(be%p%e%ip(isl),:))

       ! problem coefficients
       breg = -be%p%bs(isb)%p%s%ie(2) ! boundary region
       gamma = coeff_rob(xgs,breg)
       hb    = coeff_neu(xgs,breg)

       btypespec: select case(be%p%bs(isb)%p%btype)

        case(1) ! Neumann, Hughes flux

         ! advection coefficient
         a  = coeff_adv(xgs)
         an = matmul(transpose(a),be%p%e%n(:,isl))
         ap = 0.5_wp*(an + abs(an)) + gamma

         do l=1,base%ms
           do i=1,base%pk
             do j=1,base%pk
               aak(i,j) = aak(i,j) + wgsa(l) * &
                 base%pb(j,l,isl) * ap(l) * base%pb(i,l,isl)
             enddo
             fk(i) = fk(i) - wgsa(l) * hb(l) * base%pb(i,l,isl)
           enddo
         enddo
         
        case default
         err%lerr = .true.
         write(err%message,'(A,I5,A,I5,A,I5,A)') &
           'Unknown boundary condition type "',be%p%bs(isb)%p%btype, &
           '" on side ',be%p%bs(isb)%p%s%i,', on element ',be%p%e%i,'.'
       end select btypespec
      case default
       err%lerr = .true.
       write(err%message,'(A,I5,A,I5,A,I5,A)') &
         'Unknown boundary condition "',be%p%bs(isb)%p%bc,         &
         '" on side ',be%p%bs(isb)%p%s%i,', on element ',be%p%e%i,'.'
     end select btype
   enddo

 end subroutine cg_boundary_flux
 
!-----------------------------------------------------------------------

end module mod_cg_locmat

